####################################################################################

#Making the bourgeoisie? Values, voice, and state-provided homeownership
#Code to create Tables, Figures, and data cited in the main text and appendix 
#using Routputs created by running the previous three files (see documentation)
#NOTE: Line 122 added after original replication to correct an omission in Figure 1
#Line 99 edited to include additional specifications in Figure F4
####################################################################################

rm(list=ls())
library(AER)
library(tidyverse)
library(stringi)
library(stringr)
library(readr)
library(readxl)
library(ggplot2)
library(lubridate)
library(sandwich)
library(ri)
library(RItools)
library(xtable)
library(genderBR)
library(boot)
library(corrgram)
library(psych)
library(forcats)
library(estimatr)
library(fastDummies)
library(DeclareDesign)
library(randomizr)
library(ANOVAreplication)
library(miceadds)
library(clubSandwich)
library(texreg)
options(scipen=6)

setwd("~/Dropbox/Replication-Making-JoP/") #Set WD for replication
source("Code/functions.R") #auxiliary functions for analysis

############################################################################################
# Figure 1: ITT Effects on Values and Voice
############################################################################################

# Comparing outcomes across waves

load("Routputs/out-W1-estimates.RData")
out.all.w1 <- out.all

load("Routputs/out-W2-estimates.RData")
out.all.w2 <- out.all

#Getting ITTs attitudes
itts <- merge(out.all.w1$out.itt.lin$out.linc 
              ,out.all.w2$out.itt.lin$out.linc
              ,by=0,suffixes=c(".w1",".w2") )
rownames(itts) <- itts$Row.names
itts$Row.names <- NULL
itts <- itts[rownames(out.all.w1$out.itt.lin$out.linc),]
tp <- itts[,c(1,2,3,6,5,7,8,9,12,11,10)]

#Getting ITTs experience
exp.outcomes <- c("mobilization_index","talk_cand","client_support","turnout_2016",
                  "turnout_2018","turnout_2020","client_g",
                  "col_action_index","assoc","talk_recent",
                  "contact_gov","contact_civsoc",
                  "contact_collectively","effectiveness_index")
exp.labels <-  c("Mobilization Index","Talk to Candidate","Support Candidate"
                 ,"Turnout 2016","Turnout 2018","Exp. Turnout 2020","Clientelism"
                 ,"Col. Action Index", "Association","Talk to Politician" 
                 ,"Contact Govt.","Contact\nCiv. Soc. Org."
                 ,"Claim Collectively","Effectiveness Index")
exp.columns <- c(1,2,3,6,5,7,8,9,12,11,10)

load("Routputs/out-W2-estimates-experiences.RData")
print(comment(out.exp.all))
ittexp <- cbind(out.exp.all$out.exp.itt.lin$out.exp.lin
                ,out.exp.all$out.exp.itt.lin$out.exp.linc
                ,out.exp.all$out.exp.itt.lm$out.exp.lm
                ,out.exp.all$out.exp.itt.lm$out.exp.lmc)  
ittexp <- ittexp[exp.outcomes,]

## Perform formal comparison of estimates ##
diff.ittsc <- list()
for(i in 1:nrow(tp)){
  diff.ittsc[[i]] <- difftest(tp$Estimate.w1[i], tp$Estimate.w2[i],
                              tp$'Std. Error.w1'[i], tp$'Std. Error.w2'[i], alpha = 0.05)
}
diff.ittsc <- do.call("rbind",diff.ittsc)
rownames(diff.ittsc) <- rownames(itts)

# Voice "collective action",
voice.outcomes <- c("col_action_index","assoc","talk_recent",
                    "contact_gov","contact_civsoc",
                    "contact_collectively","effectiveness_index")
voice.labels <-  c("Voice\n(Summary)", "Association","Talk to Politician"
                   ,"Contact Govt.","Contact\nCiv. Soc. Org."
                   ,"Claim Collectively","Effectiveness\nIndex")
ittexp <- ittexp[voice.outcomes,]
plotdvoice <- ittexp[7:11]

# Definitions for Tables and Figures
main.outcomes <- c("attitudes_index","red_index","self_index","market_index")
main.labels <- c("Values\n(Summary)","Redistribution","Individualism","Market")

# Wave 1
ittw1c <- out.all.w1$out.itt.lin$out.linc[main.outcomes,c(1,2)]  

# Wave 2
ittw2c <- out.all.w2$out.itt.lin$out.linc[main.outcomes,c(1,2)]  

# Diff
ittdc <- diff.ittsc[main.outcomes,c("dom","se")]
diffsc <- cbind(ittw1c,ittw2c,ittdc)

# Report in text:
cat("The distance in ITTs between the two surveys for the values summary index is particularly striking as it varies from ",
    round(diff.ittsc[1,2],2)," to ",
    round(diff.ittsc[1,1],2)," (p-value of difference $=$",
    round(diff.ittsc[1,"pval"],3),").\n",sep="")

plotdvalues <- plotd <- diffsc #choose version with controls to plot
names(plotdvalues)[5:6] <- names(plotdvalues)[1:2] #new line (after original replic)

setEPS()
postscript("Figures/fig1.eps",
    width = 12, height = 7)
par(mfrow = c(1,2),
    oma = c(5,0,1,0) + 0.1,
    mar = c(.5,.5,.5,.5) + 0.1)
#Figure for Values with differences between waves
estcol=grep("^Estimate",names(plotdvalues))
secol=grep("^Std. Error",names(plotdvalues))
offsets <- c(.15,0,-.15) #c(0.1,-0.1)
bgs <- c(1,1,"white")
pchs <- c(24,25,21)
cols <- c(1,1,gray(0.4))
ltys <- c(1,1,2)
ys <- -1:-nrow(plotdvalues)
xs <- c(min(-1,min(plotdvalues[,estcol]-1.96*plotdvalues[,secol])),
        max(1,max(plotdvalues[,estcol]+1.96*plotdvalues[,secol])))
xs <- c(-.5,.5)
par(mar=c(0,6.5,1,1))
plot(xs,c(-0.5,-nrow(plotdvalues)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Values"
)
axis(side=2,at=ys
     ,labels=main.labels
     ,las=2,tick=F,cex.axis=1)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-7.5),lwd=2,col=gray(0.3))
for(i in 1:length(estcol)){
  segments(x0= plotdvalues[,estcol[i]]-1.96*plotdvalues[,secol[i]],
           x1= plotdvalues[,estcol[i]]+1.96*plotdvalues[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],lty=ltys[i])
  points(plotdvalues[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
estcol=grep("^Estimate",names(plotdvoice))
secol=grep("^Std. Error",names(plotdvoice))
offsets <-0
bgs <- 1
pchs <- 25
cols <- 1
ys <- -1:-nrow(plotdvoice)
xs <- c(min(-1,min(plotdvoice[,estcol]-1.96*plotdvoice[,secol])),
        max(1,max(plotdvoice[,estcol]+1.96*plotdvoice[,secol])))
xs <- c(-.5,.5)
par(mar=c(0,7,1,0.5))
plot(xs,c(-0.5,-nrow(plotdvoice)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Voice"
)
axis(side=2,at=ys
     ,labels=voice.labels
     ,las=2,tick=F,cex.axis=1)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= plotdvoice[,estcol[i]]-1.96*plotdvoice[,secol[i]],
           x1= plotdvoice[,estcol[i]]+1.96*plotdvoice[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(plotdvoice[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
#General
title(xlab = "Intent-to-treat effects",
      ylab = "",
      outer = TRUE, line = 2)
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
legend(x="bottom"
       ,legend=c("Recent lotteries survey ITT","Early lotteries survey ITT","Distance between ITTs")
       ,pch=c(24,25,21),col=cols,horiz=T
       ,pt.bg=c(1,1,"white"),bty="n",cex=1,xpd=NA,lty=c(1,1,2)
       ,inset=c(0,0))
dev.off()

############################################################################################
# Figure 2: ITT Effects on Well-Being, Expectations, and Happiness
############################################################################################

sec.outcomes <- c("welfare_index","pea","inc_above1mw","placement",
                  "prosp_self","hap_index")
sec.labels <- c("Well-Being\n(Summary)","Econ. active","Income\n>1 Min. Wage","Self placement",
                "Expectations","Happiness")

# Wave 1
ittw1c <- out.all.w1$out.itt.lin$out.linc[sec.outcomes,c(1,2)]  

# Wave 2
ittw2c <- out.all.w2$out.itt.lin$out.linc[sec.outcomes,c(1,2)]  

# Diff
ittdc <- diff.ittsc[sec.outcomes,c("dom","se")]
colnames(ittdc) <- c("Estimate Diff","Std. Error")

diffsc <- cbind(ittw1c,ittw2c,ittdc)

plotd <- diffsc
setEPS()
postscript("Figures/fig2.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(plotd))
secol=grep("^Std. Error",names(plotd))
offsets <- c(.15,0,-.15) #c(0.1,-0.1)
bgs <- c(1,1,"white")
pchs <- c(24,25,21)
cols <- c(1,1,gray(0.4))
ltys <- c(1,1,2)
ys <- -1:-nrow(plotd)
xs <- c(min(-1,min(plotd[,estcol]-1.96*plotd[,secol])),
        max(1,max(plotd[,estcol]+1.96*plotd[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,8,4,1))
plot(xs,c(-0.5,-nrow(plotd)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-treat effects"
     #,main="Early Lotteries"
)
axis(side=2,at=ys
     ,labels=sec.labels
     ,las=2,tick=F,cex.axis=1.2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-7.5),lwd=2,col=gray(0.3))
for(i in 1:length(estcol)){
  segments(x0= plotd[,estcol[i]]-1.96*plotd[,secol[i]],
           x1= plotd[,estcol[i]]+1.96*plotd[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],lty=ltys[i])
  points(plotd[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
legend(x="topright",legend=c("Recent lotteries survey ITT","Early lottery survey ITT","Distance between surveys")
       # , title = "ITT"
       ,pch=c(24,25,21),col=cols,pt.bg=c(1,1,"white"),ncol=1,bty="n"
       ,cex=1.1,xpd=NA,lty=c(1,1,2),inset=c(0,-0.1))
dev.off()


############################################################################################
# Figure 3: Heterogeneous treatment effects given wait time
############################################################################################

mod.outcomes <- c("attitudes_index","red_index","self_index","market_index",
                  "welfare_index","pea","inc_above1mw","placement",
                  "prosp_self","hap_index")
mod.labels <- c("Values\n(Summary)","Redistribution","Individualism","Market",
                "Well-being\n(Summary)","Econ. active","Income\n> 1 Min. Wage","Self placement",
                "Expectations","Happiness")


load("Routputs/out-W1-estimates-moderators.RData")
print(comment(out.mod.all))
tp.mod <- out.mod.all$out.wait[mod.outcomes,]
tp.mod$diff <- tp.mod$diff * -1 #invert so that plot reflects long wait - short wait, which is more intuitive to report

setEPS()
postscript("Figures/fig3.eps",
           width = 7, height = 10)
estcol=c(grep("^Estimate",names(tp.mod)),grep("^diff",names(tp.mod)))[c(2,1,3)]#this makes long come first
secol=c(grep("^Std. Error",names(tp.mod)),grep("^se_diff",names(tp.mod)))[c(2,1,3)]
offsets <- c(.15,0,-.15) #c(0.1,-0.1)
bgs <- c(1,1,"white")
pchs <- c(24,25,21)
cols <- c(1,1,gray(0.4))
ltys <- c(1,1,2)
ys <- -1:-nrow(tp.mod)
xs <- c(min(-1,min(tp.mod[,estcol]-1.96*tp.mod[,secol])),
        max(1,max(tp.mod[,estcol]+1.96*tp.mod[,secol])))
xs <- c(-1,1)
par(mar=c(4,8,4,1))
plot(xs,c(-0.5,-nrow(tp.mod))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-treat effects"
     #,main="Recent Lotteries"
)
axis(side=2,at=ys
     ,labels=mod.labels 
     ,las=2,tick=F,cex.axis=1.2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-4.5,-8.5,-9.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp.mod[,estcol[i]]-1.96*tp.mod[,secol[i]],
           x1= tp.mod[,estcol[i]]+1.96*tp.mod[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],lty=ltys[i])
  points(tp.mod[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
legend(x="topright",legend=c("Long wait - Control","Short wait - Control","Long wait - Short wait")
       ,pch=c(24,25,21),col=cols,pt.bg=c(1,1,"white"),ncol=1
       ,cex=1.1,xpd=NA,lty=c(1,1,2),bty="n",inset=c(0,-0.1))
dev.off()


#################### Online Appendix Analyses

#Figure A1: not created in R

#Table A1: Loading anonymized data
load("Routputs/out-a1.RData")
print(comment(out.a1))
print.xtable(xtable(out.a1$outa1a, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-a1a.tex")
print.xtable(xtable(out.a1$outa1b, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-a1b.tex")
#The remainder of Table A1 was done manually by the authors

#Table C.1
load("Routputs/out-W2-sum.RData")

sumtab <- out.sum.W2$sumtab.raw
xsumtab <- xtable(sumtab,digits=c(0,0,2,2,2,2,2,0))
label(xsumtab) <- paste("tab:summary",sep="")
caption(xsumtab) <- paste("Descriptive Statistics",sep="")
print(xsumtab,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tb-c1-raw.tex",sep=""))

sumtab <- out.sum.W2$sumtab.std
xsumtab <- xtable(sumtab, digits=c(0,0,2,2,2,2,2,0))
label(xsumtab) <- paste("tab:summary",sep="")
caption(xsumtab) <- paste("Descriptive Statistics (Standardized)",sep="")
#This table is no longer reported (raw values are more informative)
#print(xsumtab,
#      type="latex",
#      hline.after=1,
#      caption.placement="top",
#      file=paste("Tables/tb-c1-std.tex",sep=""))

#Table C.2
load("Routputs/out-W1-sum.RData")
sumtab <- out.sum.W1$sumtab.raw
xsumtab <- xtable(sumtab,digits=c(0,0,2,2,2,2,2,0))
label(xsumtab) <- paste("tab:summary",sep="")
caption(xsumtab) <- paste("Descriptive Statistics",sep="")
print(xsumtab,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tb-c2-raw.tex",sep=""))

sumtab <- out.sum.W1$sumtab.std
xsumtab <- xtable(sumtab, digits=c(0,0,2,2,2,2,2,0))
label(xsumtab) <- paste("tab:summary",sep="")
caption(xsumtab) <- paste("Descriptive Statistics (Standardized)",sep="")
#This table is no longer reported (raw values are more informative)
#print(xsumtab,
#      type="latex",
#      hline.after=1,
#      caption.placement="top",
#      file=paste("Tables/tb-c2-std.tex",sep=""))

#Table D.1

load("Routputs/out-W1-balancepooled.RData")
print.xtable(xtable(pol.bal.stats$wald, 
             hline.after = 0, only.contents = TRUE), type = "latex", 
file = "Tables/tb-d1.tex")
pol.bal.stats$p.permutation #add p-value to table

#Table D.2

load("Routputs/out-W2-balancepooled.RData")
print.xtable(xtable(t(as.matrix(pol.bal.stats$wald)), 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d2.tex")
pol.bal.stats$p.permutation #add p-value to table


#Table D.3
load("Routputs/out-W1-balancepooled.RData")
bal.labels <- c(
  "Sex (male)",
  "Race (white)", 
  "Religion (any)",
  "Children (N)",
  "Schooling (years)",
  "Registry (Cadunico)",
  "Years in Formal Employment",
  "Avg. Formal Wages")

tb3 <- bind_rows(pol.bal.stats$balage, pol.bal.stats$balnoage) %>% dplyr::select(Estimate, `Std. Error`, 
                                                                          `t value`, `Pr(>|t|)`, 
                                                                          N) 
rownames(tb3) <- c(bal.labels, paste0(bal.labels, "noage"))
print.xtable(xtable(tb3, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d3.tex")

#Table D.4 
load("Routputs/out-W2-balancepooled.RData")
controls <- pre.treat <-  c("idade", 
                            "sexor", 
                            "race_bin", 
                            "cadunico",
                            "formal_pre",
                            "av_prel")

bal.labels <- c("Age",
                "Sex (male)",
                "Race (white)", 
                "Registry (Cadunico)",
                "Years in Formal Employment", 
                "Avg. Formal Wages")

tb4 <- bind_rows(pol.bal.stats$bal, pol.bal.stats$balr) %>% dplyr::select(Estimate, `Std. Error`, 
                                                                          `t value`, `Pr(>|t|)`, 
                                                                          N) 
rownames(tb4) <- c(bal.labels, paste0(bal.labels, "femodel"))
print.xtable(xtable(tb4, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d4.tex")


#Table D.5
load("Routputs/out-W1-balancebylottery.RData")

tb5 <- bind_cols(bal.stats$edital17.2016$wald, bal.stats$edital20.2016$wald)
print.xtable(xtable(tb5, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d5.tex")
bal.stats$edital17.2016$p.permutation #add p-value to table
bal.stats$edital20.2016$p.permutation #add p-value to table


#Table D.6
load("Routputs/out-W2-balancebylottery.RData")

tb6 <- bind_cols(bal.stats$edital03.2011$wald, bal.stats$edital06.2011$wald)
print.xtable(xtable(tb6, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d6.tex")
bal.stats$edital03.2011$p.permutation #add p-value to table
bal.stats$edital06.2011$p.permutation #add p-value to table


#Table D.7
load("Routputs/out-W1-balancebylottery.RData")
tb7 <- bind_rows(bal.stats$edital17.2016$balage, bal.stats$edital17.2016$balnoage) %>% dplyr::select(Estimate, `Std. Error`, 
                                                                  `t value`, `Pr(>|t|)`, N)

bal.labels1 <- c(
  "Sex (male)",
  "Race (white)", 
  "Religion (any)",
  "Children (N)",
  "Schooling (years)",
  "Registry (Cadunico)",
  "Years in Formal Employment",
  "Avg. Formal Wages")                                                                                                                                           

rownames(tb7) <- c(bal.labels1, paste0(bal.labels1, "noage"))
print.xtable(xtable(tb7, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d7.tex")

#Table D.8

tb8 <- bind_rows(bal.stats$edital20.2016$balage, bal.stats$edital20.2016$balnoage) %>% dplyr::select(Estimate, `Std. Error`, 
                                                                                              `t value`, `Pr(>|t|)`, 
                                                                                              N) 
rownames(tb8) <- c(bal.labels1, paste0(bal.labels1, "noage"))
print.xtable(xtable(tb8, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d8.tex")

#Table D.9
load("Routputs/out-W2-balancebylottery.RData")
tb9 <- bind_rows(bal.stats$edital03.2011$balr, bal.stats$edital06.2011$balr) %>% dplyr::select(Estimate, `Std. Error`, 
                                                                                              `t value`, `Pr(>|t|)`, 
                                                                                              N) 
rownames(tb9) <- c(bal.labels, paste0(bal.labels, "0611"))
print.xtable(xtable(tb9, 
                    hline.after = 0, only.contents = TRUE), type = "latex", 
             file = "Tables/tb-d9.tex")


#Table E.1
load("Routputs/out-W1-attrition.RData")

print.xtable(xtable(attrition.stats$wald_interviewed), type = "latex", 
             file = "tables/tb-e1a.tex")
print.xtable(xtable(attrition.stats$wald_interviewed_pickedup), type = "latex", 
             file = "tables/tb-e1b.tex")

#Table E.2
load("Routputs/out-W2-attrition.RData")

print.xtable(xtable(attrition.stats.W2$wald_interviewed), type = "latex", 
             file = "tables/tb-e2a.tex")
print.xtable(xtable(attrition.stats.W2$wald_interviewed_pickedup), type = "latex", 
             file = "tables/tb-e2b.tex")


#Table E.3
texreg(list(attrition.stats$attrition_overall, attrition.stats$attrition_overall_fe, 
            attrition.stats$attrition_answered, attrition.stats$attrition_answered_fe), 
       include.ci = FALSE, 
       file =  "tables/tb-e3.tex")

#Table E.4 
texreg(list(attrition.stats.W2$attrition_overall, attrition.stats.W2$attrition_overall_fe, 
            attrition.stats.W2$attrition_answered, attrition.stats.W2$attrition_answered_fe), 
       include.ci = FALSE, 
       file =  "tables/tb-e4.tex")


#Table E.5
load("Routputs/out-W2-attrition-admin.RData")

benchmarkW2 <- texreg(list(attrition.stats.W2.admin$admin_av, 
                           attrition.stats.W2.admin$admin_formal, 
                           attrition.stats.W2.admin$survey_av, 
                           attrition.stats.W2.admin$survey_formal), include.ci = F)

print(benchmarkW2,
      type="latex",
      hline.after=1,
      caption.placement="top",
      file=paste("Tables/tb-e5.tex",sep=""))


#Data cited 
difftest(as.numeric(attrition.stats.W2.admin$survey_av$coefficients), as.numeric(attrition.stats.W2.admin$admin_av$coefficients), 
         as.numeric(attrition.stats.W2.admin$survey_av$std.error), as.numeric(attrition.stats.W2.admin$admin_av$std.error))

difftest(as.numeric(attrition.stats.W2.admin$survey_formal$coefficients), as.numeric(attrition.stats.W2.admin$admin_formal$coefficients), 
         as.numeric(attrition.stats.W2.admin$survey_formal$std.error), as.numeric(attrition.stats.W2.admin$admin_formal$std.error))


#Figure F.1
allitems <- c("tax_item","red_abstract","red_concrete",
              "effortbetter","trustothers","successalone","moneyimportant",
              "govind","competition"  )
allitems.labels <- c("Tax & Spend","Abstract Redist.","Concrete Redist.",
                     "Effort","Trust","Success","Money",
                     "Individuals v Gov.","Competition")

load("Routputs/out-W1-estimates-unplanned.RData")
out.allitems.w1 <- out.un.all$out.itt.un.lin$out.linc[allitems,]

# Wave 2
load("Routputs/out-W2-estimates-uplanned.RData")
out.allitems.w2 <- out.un.all$out.itt.un.lin$out.linc[allitems,]

setEPS()
postscript("Figures/fig-f1.eps",
           width = 14, height = 10)
par(mfrow = c(1,2),
    oma = c(5,11.5,0,0) + 0.1,
    mar = c(1,.35,1,.35) + 0.1)
#General settings
offsets <- 0#c(.22,.1,-.1,-.22) #c(0.1,-0.1)
bgs <- 1#c(1,0,1,0)
pchs <- 21#c(21,21,24,24)
cols <- 1#c(1,1,1,1)
ys <- -1:-nrow(out.allitems.w1)
#Wave 1
tp <- out.allitems.w1
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-(nrow(tp)+.25))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Recent Lotteries")
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-3.5,-7.5),col=gray(0.7))
abline(h=c(-1.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
#Wave 2
tp <- out.allitems.w2
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-(nrow(tp)+.25))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Early Lotteries")
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-3.5,-7.5),col=gray(0.7))
abline(h=c(-3.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
#general
title(xlab = "Treatment Effects",
      ylab = "",
      outer = TRUE, line = 0)
axis(side=2,at=ys
     ,labels=allitems.labels
     ,las=2,tick=F,cex.axis=1.05,outer = TRUE)
dev.off()

#F.2(a) 
itts <- cbind(out.un.all$out.itt.un.lm$out.lm 
              ,out.un.all$out.itt.un.lin$out.lin)  
itts <- itts[,c(1,2,3,6,5,7,8,9,12,11,10)] 
satisf.items <- grep("satisf",rownames(itts))
tp <- itts[satisf.items,]

the.labels <-  c("Effort Better", "Trust Others", "Success Alone", 
                 "Money Important", 
                 "Tax Item", "Redistribution (Abstract)", 
                 "Redistribution (Concrete)","Govt vs. Individual", 
                 "Comp. is Good",
                 "Expect. Self","Expect. Socio",
                 "Claim-making Govt.", 
                 "Claim-making CSO", "Effectiv. of Claim", "Effectiv. of Change",
                 "Claim - Transport.", "Claim - Health/Educ.", 
                 "Claim - Infra.", "Claim - Security",
                 "Rel. Econ. Situation", "Avg. Income", "Formal Job", 
                 "Commute (dummy)", 
                 "Help from network", "Network size (dummy)",
                 "Retrosp.","Talent Effort","Government",
                 "Satisf. Home (index)",
                 "Satisf. Size",
                 "Satisf. Conf.",
                 "Satisf. Cost",
                 "Satisf. N'hood (index)",
                 "Satisf. Violence",
                 "Satisf. Transp.",
                 "Satisf. Health",
                 "Satisf. Educ",
                 "Better Life","Would Move"
)

setEPS()
postscript("Figures/fig-f2a.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.1,-0.1)
ys <- -1:-nrow(tp)+0.5
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.4,.4)
par(mar=c(6,8,2,1))
plot(xs,c(0,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects"
     ,main="Early Lotteries")
axis(side=2,at=ys
     ,labels=the.labels[satisf.items] 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
#abline(h=c(-10.5,-12.5,-17.5,-22.5),lwd=2,col=gray(0.3))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=0+i,pch=20+i,cex=0.8,col=0+i)	
}
legend(x="bottom",legend=c("ITT (via lm)","ITT (via lin)")
       ,pch=21:22,col=1:2,pt.bg=1:2,ncol=2,bty="n",cex=0.8,inset=-0.135,xpd=NA)
dev.off()

#F.2(b)

contact.items <- grep("contact_",rownames(itts))[-c(1,2)]
tp <- itts[contact.items,]

setEPS()
postscript("Figures/fig-f2b.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
offsets <- c(0.1,-0.1)
ys <- -1:-nrow(tp)+0.5
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.4,.4)
par(mar=c(6,8,2,1))
plot(xs,c(0,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Treatment Effects"
     ,main="Early Lotteries")
axis(side=2,at=ys
     ,labels=the.labels[contact.items] 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
#abline(h=c(-10.5,-12.5,-17.5,-22.5),lwd=2,col=gray(0.3))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=0+i,pch=20+i,cex=0.8,col=0+i)	
}
legend(x="bottom",legend=c("ITT (via lm)","ITT (via lin)")
       ,pch=21:22,col=1:2,pt.bg=1:2,ncol=2,bty="n",cex=0.8,inset=-0.135,xpd=NA)
dev.off()

#Figure F.3

main.outcomes <- c("attitudes_index","red_index","self_index","market_index")
main.labels <- c("Values\n(Summary)","Redistribution","Individualism","Market")
the.columns <- c(1,2,3,6,5,7,8,9,12,11,10)

# Wave 1
ittw1 <- cbind(out.all.w1$out.itt.lin$out.lin
               ,out.all.w1$out.itt.lin$out.linc
               ,out.all.w1$out.itt.lm$out.lm
               ,out.all.w1$out.itt.lm$out.lmc)  
ittw1 <- ittw1[main.outcomes,]
tp <- ittw1

# Wave 2
ittw2 <- cbind(out.all.w2$out.itt.lin$out.lin
               ,out.all.w2$out.itt.lin$out.linc
               ,out.all.w2$out.itt.lm$out.lm
               ,out.all.w2$out.itt.lm$out.lmc)  
ittw2 <- ittw2[main.outcomes,]

setEPS()
postscript("Figures/fig-f3.eps",
           width = 14, height = 10)
par(mfrow = c(1,2),
    oma = c(5,11.5,0,0) + 0.1,
    mar = c(1,.35,1,.35) + 0.1)
#General settings
offsets <- c(.22,.1,-.1,-.22) 
bgs <- c(1,"white",1,"white")
pchs <- c(21,21,24,24)
cols <- c(1,1,1,1)
ys <- -1:-nrow(tp)
#Wave 1
tp <- ittw1
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-(nrow(tp)+.25))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Recent Lotteries")
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-1.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
#Wave 2
tp <- ittw2
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-(nrow(tp)+.25))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Early Lotteries")
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-3.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
#general
title(xlab = "Treatment Effects",
      ylab = "",
      outer = TRUE, line = 0)
axis(side=2,at=ys
     ,labels=main.labels
     ,las=2,tick=F,cex.axis=1.2,outer = TRUE)
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
legend(x=0.5,y=0.05
       ,legend=c("Lin","Lin + ctrls","Fixed effects","Fixed effects + ctrls")
       ,pch=c(21,21,24,24),col=1,pt.bg=c(1,0,1,0),ncol=2,bty="n",cex=1.2,xpd=NA)
dev.off()

#F.4
voice.outcomes <- c("col_action_index","assoc","talk_recent",
                    "contact_gov","contact_civsoc",
                    "contact_collectively","effectiveness_index")
voice.labels <-  c("Voice\n(Summary)", "Association","Talk to Politician"
                   ,"Contact Govt.","Contact\nCiv. Soc. Org."
                   ,"Claim Collectively","Effectiveness\nIndex")
ittexp <- ittexp[voice.outcomes,]

setEPS()
postscript("Figures/fig-f4.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(ittexp))
secol=grep("^Std. Error",names(ittexp))
offsets <- c(.22,.1,-.1,-.22) 
bgs <- c(1,0,1,0)
pchs <- c(21,21,24,24)
cols <- c(1,1,1,1)
ys <- -1:-nrow(ittexp)
xs <- c(min(-1,min(ittexp[,estcol]-1.96*ittexp[,secol])),
        max(1,max(ittexp[,estcol]+1.96*ittexp[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,9,2,.5))
plot(xs,c(0,-nrow(ittexp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-treat effects"
     ,main="Early lotteries survey")
axis(side=2,at=ys
     ,labels=voice.labels
     ,las=2,tick=F,cex.axis=1.2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= ittexp[,estcol[i]]-1.96*ittexp[,secol[i]],
           x1= ittexp[,estcol[i]]+1.96*ittexp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(ittexp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
legend(x="top",legend=c("Lin","Lin + ctrls","Fixed effects","Fixed effects + ctrls")
       ,pch=c(21,21,24,24),col=1,pt.bg=c(1,0,1,0),ncol=2,bty="n",cex=1.2,xpd=NA)
dev.off()

#F.5
load("Routputs/out-W2-estimates-ipw.RData")

ittsipw <- rbind(out.all$out.itt.ipw, out.all$out.itt.ipw_m)
tp <- ittsipw

the.labels.all <- c("Attitudes Index","Self-reliance",
                    "Market Beliefs","Redistribution",
                    "Welfare Index"
                    ,"Placement"
                    ,"PEA"
                    ,"Inc.>1MW"
                    ,"Happiness"
                    ,"Expectations", 
                    "Col. Action Index", "Association","Talk to Politician"
                    ,"Petition Govt.","Petition Civ. Soc. Org."
                    , "Claim Collectively","Effectiveness Index")


setEPS()
postscript("Figures/fig-f5.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-.4,.4)
par(mar=c(4,8,2,1))
plot(xs,c(0,-nrow(tp))
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="ITTs"
     ,main="IPW Analysis")
axis(side=2,at=ys
     ,labels=the.labels.all 
     ,las=2,tick=F,cex.axis=0.8)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys,lty=3,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys)
  points(tp[,estcol[i]],ys,bg=0+i,pch=20+i,cex=0.8,col=0+i)	
}
legend(x="top",legend=c("ITT+Ctrls (via IPW)")
       ,pch=21,col=1,pt.bg=1,ncol=2,bty="n",
       cex=0.8,inset=-0.01,xpd=NA)
dev.off()

# Table F.1
main.outcomes <- c("attitudes_index","red_index","self_index","market_index", 
                   "welfare_index", "placement", "pea", "inc_above1mw","hap_index","prosp_self")
w1.values <- out.all.w1$out.itt.lin$out.linc[main.outcomes,]
w2.values <- out.all.w2$out.itt.lin$out.linc[main.outcomes,]

exp.outcomes <- c("col_action_index","assoc","talk_recent",
                  "contact_gov","contact_civsoc",
                  "contact_collectively","effectiveness_index")

exp.labels <-  c("Col. Action Index", "Association","Talk to Politician"
                 ,"Contact Govt.","Contact\nCiv. Soc. Org."
                 ,"Claim Collectively","Effectiveness Index")

w2.voice <- out.exp.all$out.exp.itt.lin$out.exp.linc[exp.outcomes,]
w1.voice <- w2.voice * NA #create placeholder table 

w2.all <- rbind(w2.values,w2.voice)
w1.all <- rbind(w1.values,w1.voice)

main.labels <- c("Values\n(Summary)","Redistribution","Individualism","Market", 
                 "Well-being\n(Summary)", "Self placement", "Econ. active", "Income above 1 MW",
                 "Happiness","Expectations")

the.tab <- cbind(w1.all,NA,w2.all)[,c(1,2,3,6,4,7,8,9,10,13,11)]
rownames(the.tab) <- c(main.labels,exp.labels)
the.xtab <- xtable(the.tab,digits=3)
caption(the.xtab) <- "Point Estimates"
label(the.xtab) <- "tab:pointestimates"

print(the.xtab,
      type="latex",
      hline.after=1,
      caption.placement="top",
      digits=3,
      file=paste("Tables/tb-f1.tex",sep=""))

#Figure G.1 
main.outcomes <- c("attitudes_index","red_index","self_index","market_index")
the.labels <- c("Values\n(Summary)","Redistribution","Individualism","Market")
the.columns <- c(1,2,3,6,5,7,8,9,12,11,10)

# Wave 1
cacew1 <- cbind(out.all.w1$out.cace.lm$out.lmcace
                ,out.all.w1$out.cace.lm$out.lmcacec)  
cacew1 <- cacew1[main.outcomes,the.columns]


# Wave 2
cacew2 <- cbind(out.all.w2$out.cace.lm$out.lmcace
                ,out.all.w2$out.cace.lm$out.lmcacec)  
cacew2 <- cacew2[main.outcomes,the.columns]


setEPS()
postscript("Figures/fig-g1.eps",
           width = 14, height = 10)
par(mfrow = c(1,2),
    oma = c(5,11.5,0,0) + 0.1,
    mar = c(1,.35,1,.35) + 0.1)
#General settings
bgs <- c(1,0)
pchs <- c(24,25)
offsets <- c(0.1,-0.1)
#Wave 1
tp <- cacew1
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-1,1)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-nrow(tp)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Recent lotteries")
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=-4.5,col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)	
}
#Wave 2
tp <- cacew2
estcol=grep("^Estimate",names(tp))
secol=grep("^Std. Error",names(tp))
ys <- -1:-nrow(tp)
xs <- c(min(-1,min(tp[,estcol]-1.96*tp[,secol])),
        max(1,max(tp[,estcol]+1.96*tp[,secol])))
xs <- c(-1,1)
par(mar=c(4,2,2,1))
plot(xs,c(-0.5,-nrow(tp)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab=""
     ,main="Early lotteries"
)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-4.5,-8.5,-9.5),col=gray(0.7))
for(i in 1:length(estcol)){
  segments(x0= tp[,estcol[i]]-1.96*tp[,secol[i]],
           x1= tp[,estcol[i]]+1.96*tp[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i])
  points(tp[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=1)		
}
#general
title(xlab = "Treatment effects on compliers",
      ylab = "",
      outer = TRUE, line = 0)
axis(side=2,at=ys
     ,labels=the.labels
     ,las=2,tick=F,cex.axis=1.2,outer = TRUE)
par(mfrow=c(1, 1), oma=rep(0, 4), mar=rep(0, 4), new=TRUE)
plot(0:1, 0:1, type="n", xlab="", ylab="", axes=FALSE)
legend(x=0.5,y=0#x="bottom"
       ,legend=c("CACE","CACE+Cltrs")
       ,pch=24:25,col=1,pt.bg=c(1,0),ncol=2,bty="n",cex=1.2,inset=-0.005,xpd=NA)
dev.off()

#Appendix H data cited
load("Routputs/allw.Rda")

part <- subset(allw,select=c(wave2,edital03.2011,edital06.2011,edital17.2016,edital20.2016))
part$editais2011 <- part$edital03.2011 | part$edital06.2011
part$editais2016 <- part$edital17.2016 | part$edital20.2016

table(Wave2=part$wave2,Editais.2011=part$editais2011)
table(Wave2=part$wave2,Editais.2016=part$editais2016)
sum(part$wave2==1&part$editais2016==T)

#Comparing samples
cat("Of the ",sum(part$wave2)," observations in the early lotteries survey ",
    sum(part$wave2==1&part$editais2016==T)," (",
    round(100*sum(part$wave2==1&part$editais2016==T)/sum(part$wave2),1),"%) correspond to individuals who also participated in the late lotteries. Of the ",
    sum(part$wave2==0)," observations in the late lotteries survey ",
    sum(part$wave2==0&part$editais2011==T)," (",
    round(100*sum(part$wave2==0&part$editais2011==T)/sum(part$wave2==0),1),"%) correspond to individuals who also participated in the early lotteries.\n",sep="")

load("Routputs/out-inscritosearlylate.RData")
print(comment(insc))
sumtab <- round(100*prop.table(table(l2011=insc$w2,l2016=insc$w1,useNA="always")),1)

cat("A total of "
    ,nrow(insc)," individuals participated in at least one of the two pairs of lotteries from which we drew our sample."
    ," Of these, "
    ,sumtab[1,2],"% participated only in the 2011 lotteries, "
    ,sumtab[2,1],"% participated only in the 2016 lotteries, and "
    ,sumtab[1,1],"% participated in both."
    ," This means that "
    ,round(prop.table(table(subset(insc,w2==T)$w1,useNA="always"))[1]*100,1),"% of those in the 2011 lotteries also participated in the 2016 lotteries while "
    ,round(prop.table(table(subset(insc,w1==T)$w2,useNA="always"))[1]*100,1),"% of those in the 2016 lotteries were also in the earlier ones.",sep="")


#F-test
load("Routputs/out-bal-waves.RData")
out.bal.waves$ftest

#Figure H.1 
bal.labels <- c("Age",
                "Sex (male)",
                "Race (white)",
                "Registry (Cadunico)",
                "Years in Formal \n Employment", 
                "Avg. Formal Wages")

ballmr <- out.bal.waves$bal

setEPS()
postscript("Figures/fig-h1.eps",
           width = 7, height = 7)
par(mar=c(4,8,1.5,1))
plot(ballmr[,"Std. Estimate"],-1*1:nrow(ballmr)
     ,bty="n",xlab="Standardized Mean Differences",ylab="",
     xlim=c(-.5,.5),yaxt="n",
     pch=ifelse(ballmr[,4]>0.1,22,24),xpd=NA,
     main=paste("Balance Between Waves"))
abline(v=c(-.25,0,.25),lty=c(2,1,2),col=c(gray(0.5),1,gray(0.5)))
axis(side=2,at=-1*1:nrow(ballmr),labels=bal.labels
     ,tick=F,las=2,line=-2,xpd=NA,cex.axis=0.9)
legend(x="topright",legend=c("p > 0.1","p < 0.1"),pch=c(22,24)
       ,bg=c(0,gray(0.5)),horiz=F,inset=0.05,xpd=NA)
dev.off()

#Figure H.2

ballmro <- out.bal.waves$balo

outcomes.labels <-  c("Effort Better", "Trust Others", "Success Alone", 
                      "Money Important", 
                      "Redistribution (Abstract)", 
                      "Redistribution (Concrete)","Govt vs. Individual", 
                      "Comp. is Good"
                      ,"PEA"
                      ,"Inc.>1MW"
                      ,"Self-placement"
                      ,"Pos. Affect"
                      ,"Blue Affect"
                      ,"Cantril"
                      ,"Expectations")

setEPS()
postscript("Figures/fig-h2.eps",
           width = 7, height = 7)
par(mar=c(4,8,1.5,1))
plot(ballmro[,"Std. Estimate"],-1*1:nrow(ballmro)
     ,bty="n",xlab="Standardized Mean Differences",ylab="",
     xlim=c(-.5,.5),yaxt="n",
     pch=ifelse(ballmro[,4]>0.1,22,24),xpd=NA,
     main=paste("Control Group Raw Outcomes Between Waves"))
abline(v=c(-.25,0,.25),lty=c(2,1,2),col=c(gray(0.5),1,gray(0.5)))
axis(side=2,at=-1*1:nrow(ballmro),labels=outcomes.labels
     ,tick=F,las=2,line=-2,xpd=NA,cex.axis=0.9)
legend(x="topright",legend=c("p > 0.1","p < 0.1"),pch=c(22,24)
       ,bg=c(0,gray(0.5)),horiz=F,inset=0.05,xpd=NA)
dev.off()

#Table H.1 
load("Routputs/preds-comparing.Rda")

outcomes <- c(
  "attitudes_index",
  "red_index",
  "self_index",
  "market_index",
  "welfare_index",
  "pea",
  "inc_above1mw",
  "placement",
  "prosp_self",
  "hap_index"
)
rownames(preds) <- outcomes
xpreds <- xtable(preds)
caption(xpreds) <- "Predicting Control Group Outcomes with Pre-Treatment Covariates"
label(xpreds) <- "tab-predoutcomes"
print(xpreds,type="latex",
      hline.after=1,
      caption.placement="top",
      digits=3,
      file=paste("Tables/tb-h1.tex",sep=""))


#Figure H.3 

# Wave 1 
ittw1c <- out.all.w1$out.itt.lin$out.linc[main.outcomes,c(1,2)]  
ittw1 <- out.all.w1$out.itt.lin$out.lin[main.outcomes,c(1,2)]  

# Wave 2
load("Routputs/out-W2-matchedestimates.RData")
out.all.w2w <- out.allw #use matched (note w for weighted)
ittw2c <- out.all.w2w$out.itt.linw$out.lincw[main.outcomes,c(1,2)]  

## Perform formal comparison of estimates ##
diff.ittsc <- list()
for(i in 1:nrow(ittw1c)){
  diff.ittsc[[i]] <- difftest(ittw1c$Estimat[i], ittw2c$Estimate[i],
                              ittw1c$'Std. Error'[i], ittw2c$'Std. Error'[i], alpha = 0.05)
}
diff.ittsc <- do.call("rbind",diff.ittsc)
rownames(diff.ittsc) <- rownames(ittw2c)

# Diff
ittdc <- diff.ittsc[main.outcomes,c("dom","se")]
colnames(ittdc) <- c("Estimate Diff","Std. Error")
diffsc <- cbind(ittw1c,ittw2c,ittdc)
plotd <- diffsc 

main.labels <- c("Values\n(Summary)", "Redistribution",  "Individualism", "Market")

setEPS()
postscript("Figures/fig-h3.eps",
           width = 7, height = 10)
estcol=grep("^Estimate",names(plotd))
secol=grep("^Std. Error",names(plotd))
offsets <- c(.15,0,-.15) #c(0.1,-0.1)
bgs <- c(1,1,"white")
pchs <- c(24,25,21)
cols <- c(1,1,gray(0.4))
ltys <- c(1,1,2)
ys <- -1:-nrow(plotd)
xs <- c(min(-1,min(plotd[,estcol]-1.96*plotd[,secol])),
        max(1,max(plotd[,estcol]+1.96*plotd[,secol])))
xs <- c(-.5,.5)
par(mar=c(4,8,4,1))
plot(xs,c(-0.5,-nrow(plotd)-0.25)
     ,type="n",bty="n",yaxt="n"
     ,ylab="",xlab="Intent-to-treat effects"
)
axis(side=2,at=ys
     ,labels=main.labels
     ,las=2,tick=F,cex.axis=1.2)
abline(v=0,lty=1,col=gray(0.7))
abline(h=ys+0.5,lty=3,col=gray(0.7))
abline(h=c(-7.5),lwd=2,col=gray(0.3))
for(i in 1:length(estcol)){
  segments(x0= plotd[,estcol[i]]-1.96*plotd[,secol[i]],
           x1= plotd[,estcol[i]]+1.96*plotd[,secol[i]],
           y0=ys+offsets[i],
           y1=ys+offsets[i],lty=ltys[i])
  points(plotd[,estcol[i]],ys+offsets[i],bg=bgs[i],pch=pchs[i],cex=0.8,col=cols[i])	
}
legend(x="topright",legend=c("Recent lotteries survey ITT","Early lotteries survey ITT","Distance between surveys")
       ,pch=c(24,25,21),col=cols
       ,pt.bg=c(1,1,"white"),ncol=1,bty="n",cex=1.1,xpd=NA,lty=c(1,1,2),inset=c(0,-0.1))
dev.off()

# Figure I.1
load("Routputs/out-W1-bal-moderators.RData")
ballmr <- out.mod.ballmr$out.bal.wait
bal.labels <- c("Sex (male)", "Race (white)", "Religion (any)", "Children (N)", "Schooling (years)", "Registry (Cadunico)",          
                "Years in Formal \n Employment", "Avg. Formal Wages")  

setEPS()
postscript("Figures/fig-i1.eps",
           width = 7, height = 7)
par(mar=c(4,8,1.5,1))
plot(ballmr[,"Std. Estimate"],-1*1:nrow(ballmr)
     ,bty="n",xlab="Standardized Mean Differences",ylab="",
     xlim=c(-.5,.5),yaxt="n",
     pch=ifelse(ballmr[,4]>0.1,22,24),xpd=NA,
     main=paste("Balance Between Long- and Short-Wait"))
abline(v=c(-.25,0,.25),lty=c(2,1,2),col=c(gray(0.5),1,gray(0.5)))
axis(side=2,at=-1*1:nrow(ballmr),labels=bal.labels
     ,tick=F,las=2,line=-2,xpd=NA,cex.axis=0.9)
legend(x="topright",legend=c("p > 0.1","p < 0.1"),pch=c(22,24)
       ,bg=c(0,gray(0.5)),horiz=F,inset=c(-0.025,0.05),xpd=NA)
dev.off()

#Figure J.1
load("Routputs/out-cesop.RData")

all.regs <- out.all.cesop$all.regs
ev.mean <- out.all.cesop$ev.mean

setEPS()
postscript("Figures/fig-j1.eps",
           width = 12, height = 6)
par(mar=c(4,4,2,1))
the.plot(all.regs[order(ev.mean,decreasing = T)],ymin=4,ymax=8)
mtext(side=2,"Evaluation of Program",line=2.5)
legend(x="bottom",
       legend=c("0 to 1 MW","1 to 2 MW","2 to 5 MW","5 to 10 MW","> 10 MW","Average"),
       pch=c(21:25,NA),col=1,pt.bg=c(1,"white","white","white",1,NA),
       bty="n",lty=c(rep(1,5),2),horiz=F,ncol=2
       ,cex=1,inset=-0.175,xpd=NA
) 
axis(side=3,at=c(1:8),tick=F,line=-3
     ,labels=c("MCMV","Bolsa\nFamilia","PROUNI","PRONATEC","ENEM","Farmacia\nPopular","FIES","Mais\nMédicos")[order(ev.mean,decreasing = T)]
     ,cex.axis=1)
dev.off()

